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We use the Bethe Ansatz to derive analytical expressions for the current statistics in the asym- 
metric exclusion process with both forward and backward jumps. The Bethe equations are highly 
coupled and this fact has impeded their use to derive exact results for finite systems. We overcome 
this technical difficulty by a reformulation of the Bethe equations into a one variable polynomial 
problem, akin to the functional Bethe Ansatz. The perturbative solution of this equation leads to 
the cumulants of the current. We calculate here the first two orders and derive exact formulae for 
the mean value of the current and its fluctuations. 
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I. INTRODUCTION 

^ . 

["t I , The asymmetric exclusion process (ASEP) plays a seminal role in non-equilibrium physics of low dimensional 
systems pj]. In its simplest version, the ASEP describes a system of particles, randomly hopping on a lattice with 
hard-core exclusion interaction so that a lattice site can be occupied by only one particle at a given time. Due to 
I— 1- its minimal character, this model appears as a building block in many seemingly unrelated fields 2]. By virtue of 
' different mappings, the ASEP can be interpreted as a model for RNA transcription |3j], hopping conductivity, polymers 
in random media, surface growth traffic fiow, molecular motors Q etc... In the one-dimensional case, many exact 
results have been derived for the ASEP (for a review, see e.g., @, 01)- As a result, the relations between the intrinsic 
stochasticity of the dynamics, the external drive and the particles interactions are better understood. The fact that 
ASEP in one-dimension is an exactly solvable model should not be considered as just an elegant mathematical anomaly 
at odds with physical relevance. Indeed, many of the exact results obtained for ASEP have shed light on the behaviour 
■ of general driven diffusive systems by providing us with effective phenomenological descriptions that can be applied 
' to more realistic models [s^. Examples of such descriptions that stem from mathematical results are: shock fronts to 

a model boundary induced phase transitions, the interpretation of shocks as real space condensation (related to zero 
I ' range processes) and the additivity principle [ll|. Besides, the ASEP is a good toy- model to test the validity 

. of general claims about non-equilibrium systems: for example, the Gallavotti-Cohen fluctuation theorem is satisfied 
P5 by the ASEP (and by more general Markovian systems) as can be shown by elementary methods [l^] whereas the 
O , proof for deterministic dynamical systems requires some restrictive hypothesis and is far more technical. 
, Exact solutions for the exclusion process have been obtained by using several different approaches, and in particular 

the Matrix Product representation and the Bethe Ansatz. The Matrix Product representation was first introduced in 
P\| . [l^ to study the stationary state and the phase diagram of the ASEP with open boundaries. The main idea consists 
J> ' in representing the stationary state as a trace over a suitable, usually quadratic, algebra; this technique has been 
On , generalized to many different models, including systems with shock profiles and with different classes of particles 
' [fl [Tsl . [TgI . [lH . An exhaustive and pedagogical review on the Matrix method can be found in [l^ . The Bethe Ansatz 
was first used to calculate the spectral gap of the ASEP and the associated dynamical exponent [l^, [2l|, [l^, [H, [13] . 
Indeed, the Markov matrix of the ASEP that encodes its stochastic evolution can be mapped exactly onto a non- 
hermitian spin chain hamiltonian which is integrable. The Bethe Ansatz also allows to study spectral degeneracies 
[2^ , and to investigate variants of the ASEP and more general particle hopping processes 2lL, i28,j 2S\ (for a review 
see e.g., [soj). 

A particularly important physical quantity in the ASEP is the statistics of the current in the stationary regime. 
This current becomes a local height variable when the ASEP is translated into a Random Solid On Solid model that 
describes the growth of a random interface. Indeed, in this mapping, a forward random jump of a particle through 
5J] ' a bond corresponds to a random deposition event of a unit 'brick' on the interface; a backward jump corresponds to 
\ the evaporation of a brick. The time integrated current through a bond of the ASEP is therefore equivalent to the 
total height of the interface at a given point. In the continuous limit, the motion of this interface is described by 
the Kardar-Parisi-Zhang (KPZ) equation (see e.g., [3]). The exclusion process in one dimension is thus a discretized 
version of the KPZ equation and exact results about the ASEP have therefore interesting interpretations in terms of 
surface growth. 

For the exclusion process on a periodic ring, the mean value of the stationary current through a bond can be easily 
derived from elementary combinatorics; in the limit of a large system the mean current is given by the density of 
vacancies multiplied by the asymmetry rate. However, the higher moments of the current in the stationary state are 
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much more difficult to calculate. In fact, the full statistics of the current was determined only for the particular case 
of the totally asymmetric exclusion process (TASEP), where the particles are allowed to jump only in one direction. 
For any system size, an analytical expression for the cumulant generating function was obtained, leading to an exact 
formula for the large deviation function [3l|, [s^]. This result was derived using the Bethe equations which, for the 
TASEP, can be solved explicitely thanks to a decoupling property that reduces them to a one variable polynomial 
equation plus a self-consistency condition [HI, 113, HOl ■ 

In the general case, when jumps on both directions are allowed, the Bethe equations do not decouple and it has 
not been possible to use them to derive exact results for finite systems. An exact formula for the fluctuation of the 
current (i.e., the second moment of the current) in the long time limit could however be derived using an extension 
of the matrix method [H,[33]. But higher moments appeared to be out of reach. 

The aim of the present work is to derive analytical results for the current statistics in ASEP with forward and 
backward jumps (sometimes called the partially asymmetric exclusion process) from the Bethe Ansatz. We overcome 
the technical difficulty that hindered the solution of the Bethe equations in the general case by reducing them to an 
effective one variable problem thanks to a suitable reformulation, akin to the so-called functional Bethe Ansatz. This 
one variable equation can be interpreted as a purely algebraic question involving a divisibility condition between two 
polynomials. In this work, we use this formalism to derive the expressions of the mean value of the current and its 
variance. Our technique can be used to calculate the current cumulant to any desired order. 

The outline of this work is as follows. In section Ull we explain that the cumulant generating function can be 
expressed as the maximal eigenvalue of a suitable deformation of the Markov matrix where the deformation parameter 
represents the fugacity of the jumps. In section [Till we give the Bethe equations that allow to diagonalize this matrix. 
The reformulation of the Bethe equations as a problem in polynomial divisibility is done in section HVl In section IVl 
we solve perturbatively this purely algebraic problem to the second order with respect to the jump fugacity. This 
allows us to derive the exact formulae for the mean value and the variance of the current in section IVII The last 
section is devoted to concluding remarks. Some technical derivations are given in the appendices. 

II. CURRENT STATISTICS AS AN EIGENVALUE PROBLEM 

A. The asymmetric exclusion process 

The exclusion process on a periodic one dimensional lattice with L sites (sites i and L-\-i are identical) is a stochastic 
interacting particle model in which each lattice site is occupied by at most one particle at a given time {exclusion 
rule). The system evolves with time according to a stochastic dynamics: a particle on a site i at time t jumps, in 
the interval between t and t -\- dt, with probability p dt to the neighbouring site i + 1 if this site is empty and with 
probability q dt to the site i — 1 if this site is empty. The jump rates p and q are normalized such that p -\- q = 1. 
The special case where the jumps are totally biased in one direction (p = 1 or g = 1) is called the totally asymmetric 
exclusion process (TASEP). For p = q = 1/2, the exclusion process is symmetric (SEP). If the number of particles in 

the ring is N, the total number of configurations is given by the binomial coefficient 

We call Pt{C) the probability that the system is in the configuration C at time t. As the exclusion process is a 
continuous- time Markov process, the time evolution of -Pt(C) is determined by the master equation 

|Pi(C) = ^A/(C,C')Pt(C') = {Mo{C,C') + AhiCC) + M.,{C,C'))Pt{C') . (1) 
c c 

The Markov matrix M encodes the dynamics of the exclusion process: the non-diagonal element Mi{C,C') represents 
the transition rate from configuration C to C where a particle hops in the forward [i.e., anti-clockwise) direction, the 
non-diagonal element M_i(C,C') represents the transition rate from configuration C to C where a particle hops in 
the backward {i.e., clockwise) direction. The diagonal term Mo(C,C) ~ — J2c'^c {Mi{C\C) + i\/_i(C',C)) represents 
the exit rate from the configuration C. 

B. Generalized master equation for current statistics 

We call Yt the total distance covered by all the particles between time and time t and Pt{C, Y) the joint probability 
of being in the configuration C at time t with Yt = Y. An evolution equation, analogous to equation ([T]), can be 




written for Pt{C, Y) as foUows : 

j^PtiC,Y) = (^Mo{C,C')Pt{C',Y)+Ah(C,C')Pt{C',Y-l) + M_,{C,C')Pt{C',Y + l)'^ . (2) 
c 

We now recall how the full statistics of Yt can be determined [l^, [3l[. In terms of the generating function Ft{C) 
defined as 

+00 

Ft{C)= Y ^'"^PtiC,y), (3) 

F=-oo 

equation ^ takes the simpler form : 

^F,{C) = Y {Mo{C,C') + eUhiCC) + e-^M^,{C,C'))Ft{C') = ^ Af (7)(C,C')Ft(C') . (4) 
c c 

This equation is similar to the original Markov equation ([T]) for the probability distribution Pt{C) but where the 
original Markov matrix M is deformed into M{-j) which is given by 

M(7) = Mo + c'Mi + e~^M_i . (5) 

We emphasize that M{'^), that governs the evolution of Ft{C), is not a Markov matrix for 7 7^ (the sum of the 
elements in a given column does not vanish). 



C. Long time limit and maximal eigenvalue 

In the long time limit, t 00, the behaviour of Ft{C) is dominated by the largest eigenvalue A(7) of the matrix 
M(7) : 

Fi(C)^e^— (^)*(C|£;„,ax(7)), (6) 
where the ket |£'niax(7)) is the eigenvector corresponding to the largest eigenvalue. Therefore, when t —^ 00, we obtain 

( eT^' ) = ^ Ft(C) ~ e^— . (7) 
c 

More precisely, we have 

hm ilog(e^^*) =£^max(7)- (8) 

t^oo t 

The function i?rnax(7) contains the complete information about the cumulants of the total current Yt in the long time 
limit. For example, the total current J and the diffusion constant A are given by 

(Yt) di;„,ax(7) 



J = lim 



A = lim 



{Y,^)-{Ytf d^E^Ul) 



(9) 

7=0 

(10) 

7=0 



t^oo t d7^ 

Thus, the cumulants of Yt can be determined by carrying out a perturbative expansion of £^max(7) with respect to 7 
(a similar method has been used, in a different context, in [ssl]). The importance of the maximal eigenvalue i?max(7) 
of the matrix Af(7) also stems from the fact that it is closely related to the large deviation function G for the total 
current. We recall that the large deviation function G is defined as 



1 



Yt 



G(j)= lim - log Prob ^=j . (11) 

From equations ^ and ([TT]) . we find 

{e^Yt ) ^ ~ / e*(^(^)+^^)dj , (12) 



and deduce by the saddle-point approximation that the maximal eigenvalue -Emax(7) is the Legendre transform of the 
large deviation function G{j) 

E^Ul) = max (G{j) + 7j) . (13) 



D. Restatement of the problem 



We want to study the statistical properties of the total current in the partially asymmetric exclusion process. We 
have seen that in the long time limit, the maximal eigenvalue £'inax(7) of the deformed matrix Af (7) is the generating 
function of the cumulants of the current, i.e., the power-series expansion of -E'max(7) in the vicinity of 7 = allows 
us to determine the statistical properties of the current. In the following sections, we shall first explain how to 
diagonalize the Matrix Af(7) using the Bethe Ansatz; this method will allow us to write any eigenvalue of M(7) as a 
symmetric function of the roots of a system of coupled polynomial equations (the Bethe equations). Then, we shall 
develop a perturbative scheme to expand the maximal eigenvalue -E'max(7) in powers of 7, when 7 — > 0. The first 
order expansion will give us the the current J and the second order term will lead to the diffusion constant A. 



III. THE BETHE EQUATIONS 



The deformed matrix M(7) can be diagonalized by Bethe Ansatz. A vector P over the configuration space is an 
eigenvector of M (7) if it satisfies 

M{-f)P ^ E{-f)P . (14) 

By representing a configuration by the positions of the N particles on the ring, (ri, r2, . . . , r^r) with 1 < ri < r2 < 
• • • < rjv < i, the eigenvalue equation (fTi|) becomes 

EP{ri,...,rN) = 

^p[e^P(ri,...,ri_i, n - 1, n+i, . . . , tat) - P(ri, . . . , r„)] + 

i 

^q[e~''P{ri,...,rj_i, + 1, r^+i, . . . , tat) - F(ri, . . . , rjv)] , (15) 
j 

where the sum runs over the indices i such that r^-i < — 1 and over the indices j such that r j + 1 < rj^i ; these 
conditions ensure that the corresponding jumps are allowed. Following the coordinate Bethe Ansatz, we assume that 
the eigenvector P can be written in the form 

P(ri, . . . , r„) = ^ A, zl\^^ z;^(2) . . . z^^^^ , (16) 

where S„ is the group of the n\ permutations of n indices. The coefficients {A„} are rational functions of the fugacities 
{zi, . . . , z„}. The expression (jl6p represents an eigenvector of ^^(7) if {zi, . . . , z„} satisfy the Bethe equations [2ll.[30j: 

N 

L 



(_,)A.-i-rT go-^^.^.-b + g)^.+Pe^ for ^^l...N, {17} 
fJ-^qe-iz,Zj-{p + q)zj+pe'r ' 



and the corresponding eigenvalue of Af (7) is given by 

N , N 



E{r, zi,Z2 ...ZN)=pe^J2-+ 96"^ + • (1^) 



.2. 
i— 1 t—i 



For 7 = 0, we know that the maximal eigenvalue of the Markov matrix M is equal to and corresponds to the 
degenerate solution Zi — I for all i. 

Remark: The Gallavotti- Cohen Invariance. The Bethe equations ()17p and equation ()18|) are invariant under the 
transformation z ^ ^ , 7 ^ log ^ — 1- This symmetry implies that the spectrum of M(7) and that of Af(log ^ — 7) 
are identical. This functional identity is satisfied in particular by the largest eigenvalue of M and we have -Emax(7) = 
£^max(log p — 7) • This identity implies, using equation (fT3|) . that the large deviation function satifies the symmetry 



G(j) = G(-j)-(log^)j. (19) 



This relation is a special case of the general Fluctuation Theorem valid for a large class of systems far from equilibrium. 
It was derived for more general Markovian systems in [12]. 



A. A useful change of variables 



We introduce N auxiliary variables (j/i, at) defined as 



y. = l^^^ for* = 1...7V, (20) 

1 — a;e ^z,- 



where we have introduced the asymmetry parameter x: 



x=-. (21) 
P 



We remark that the change of variables (j20p is ill-defined for x = \ which corresponds to the symmetric exclusion 
process. In the following, our calculations will always be performed for x < \. Our results will extend to the symmetric 
case by taking the limit a; — ^ 1 in the final expressions. The Bethe equations (|17p now become 



L N 



^-xyij j-J-^ xy, - yj 



11 ^' for i = l...N. (22) 



These equations are simpler than the original ones because they involve only linear polynomials in the yiS. By taking 
the product of the Bethe equations ([H]) over all the values of j, we obtain 



(N \^ N 

l\l-xy,J il^xy.-y, 

This relation stems from the translation invariance of the model (momentum conservation). 
In terms of the ?/i's, the eigenvalue (fT8| reads 

When 7^0, all the roots 2/4(7) that correspond to the maximal eigenvalue -E'max(7) of M(7) converge to the 
degenerate solution lim^^o2/i = and the maximal eigenvalue of ^(7) also converges to 0. Using equation (|23p. we 
therefore find that, for small enough values of 7, the roots yi{j) satisfy the relation 



N _ 

e^^n^^^ = l- (25) 
I- xyi 



This relation, which is a simple consequence of the Bethe equations, will be useful in the following to select the Bethe 
roots that correspond to -Emax(7)- 

B. The TASEP case 

The Bethe equations (j22l) are a coupled non-linear system of polynomial equations in the variables j/i, . . . , j/at. 
Deriving exact results from these equations is a daunting task. However, for the special case of the totally asymmetric 
exclusion process (TASEP), which corresponds to p = 1 and g = a: = 0, the Bethe equations can be reduced to an 
effective one variable problem. Indeed, for a; = 0, the equations read 

e^'^(l-y.)'' = (-1)^-^=1^ • (26) 

Thus, all the Bethe roots yi are solutions of the one-variable polynomial equation 

e^'^(l-T)^ + Cr^ =0, (27) 



where the constant C must be determined self-consistently by the relation 



1 

C={-irl[-- (28) 

j=i 

This crucial 'decoupling' property of the the Bethe equations for a; = 0, has lead to an exact calculation of the TASEP 
spectral gap [lOj HH, Hj] and has allowed Derrida and Lebowitz to calculate the complete large deviation function of 
the current for any finite values of L and N [3l|. This effective decoupling also explains the spectral degeneracies 
of the TASEP Markov matrix [25j. Hence, the use of the Bethe Ansatz has been restricted mostly to TASEP (for a 
review see e.g. [SO]). 

For the partially asymmetric exclusion process, the Bethe equations are highly coupled to one another and can not 
be simply reduced to an effective one variable equation. Because of this technical difhculty for < a; < 1, it has not 
been possible to extract from the Bethe Ansatz any exact solution for finite systems. However, when L oo, the 
Bethe equations reduce to an integro-differential equation for the density of roots, which was analyzed by Kim et 
al. [H, [isl to derive the spectral gap and the current large deviation function. 



IV. REFORMULATION OF THE BETHE EQUATIONS 

We note that in the N Bethe equations (f22| all the variables yi play a similar role. This remark suggests that we 
should introduce an auxiliary variable T that plays a symmetric role with respect to all the j/i's. We suppose that T 
satisfies the following equation 



e^^ 



L N 



1 — xT J xT — y.j 



j-j^x^ for ^ = 1...N, (29) 



where the y^'s are now interpreted as parameters of the problem. This expression can be rewritten as a one variable 
polynomial equation for the unknown T: 

N N 

P{T) = with P{T) = e^T(l - Tf \{{xT - y,) + (1 - xT)"^ J|(T - xy,) . (30) 

The N Bethe equations (j22p imply that yi is a root of P{T) ioi i — 1, . . . ,N. Thanks to the auxiliary variable T, the 
Bethe equations have been reduced to an effective one variable problem with N parameters. We can now proceed as 
follows: (i) Find the roots of the polynomial P{T) with the unknown T and with N parameters yi, . . . , y^- (ii) Select 
N roots, amongst the L + N solutions of P{T) = 0, and identify these selected roots to the y^'s. This identification 
leads to N self-consistent equations (recall that for TASEP we had only one self-consistency condition). 

It is possible to perform these steps using contour integration in the complex plane as in the TASEP case [U [27l . 
[3ll.[3^. However, the calculations will be greatly simplified if the problem is formulated in a purely algebraic manner, 
as follows. Let us define the polynomial Q{T) as 

N 

QiT) = l[{T~yj). (31) 

The roots of Q are exactly the Bethe roots j/i , . . . , y^ (equivalently, Q is the generating function of the symmetric 
polynomials in j/i, . . . , y^)- The polynomial P{T), defined in equation ((50)) . can then be written as follows 

P{T) = e^^{l - TfQ{xT) + (1 - xPf x^Q . (32) 

The fact that the Bethe roots yi, . . . ,yN are roots of the polynomial P{T) implies that Q{T) divides P{T). Therefore 
there exists a polynomial R{T) of degree L such that P{T) = Q{T)R{T), i.e., such that 

Q{T)R{T) = e^^(l - TfQixT) + (1 - xTfx^Q ('^^ . (33) 



Substituting T = in this equation and taking into account that yt is a root of Q{T) we obtain 



l-xyij Q{xyi) 



(34) 



Using the expression ||HT|) for Q{T), we find that this relation is identical to the Bethe equation (|^^ . We remark that 
this reformulation of the Bethe equations as a problem of polynomial divisibility has been used in various contexts 
[3^ [37L [ssj and is closely related to the functional Bethe Ansatz [s^ [H, [s^ . 



A. Expression of the eigenvalue 

The eigenvalue i?(7), defined in equation ((24|) . can be expressed in terms of the polynomial Q{T): 



Q(l) X Q{l/x) 

This formula can be simplified with the help of the 'Q-R equation' 155]) as follows. Substituting T = 1 in the 
equation (|33p we find 

Q{l)R{l)^{l~xfx''Q{^ . (36) 
If we differentiate equation ([55)1 with respect to T and then substitute T = 1 we obtain 



Q'(l)i?(l) + Q(l)i?'(l) = -Lx"+i(l - x)^-iQ + x^-i(l - x)^g' Q 

Taking the ratio of the last two equations, we find that £^(7) can be rewritten as 

E(-i) _ Lx R'{1) 
p(l - a;) ^ ^ l-x ^ i?(l) ■ 



(37) 



(38) 



This is the expression of £'(7) that will be used in the sequel. 

Equation (|25p . that allows to select the roots yi corresponding to the maximal eigenvalue, is similarly rewritten in 
terms of Q and R as follows 



e^7^(i) ^ 1 . (39) 



Using equation p6p. an alternative form is obtained 

i?(l) =e^^(l-a;)^. (40) 
This relation will be very useful in the sequel to simplify some calculations. 

V. PERTURBATIVE SOLUTION OF THE FUNCTIONAL BETHE ANSATZ EQUATIONS 

In this section, we explain how to solve equation (|33p order by order in 7 for the roots yi that correspond to the 
maximal eigenvalue of the matrix M(7). 

We first develop the polynomials Q and R in powers of 7 

N 00 

Q{T) = \{{T~ yj) = I'^QniT) = Oo(T) + 7Qi(T) + i^Q2{T) + . . . (41) 

CSO 

R{T) = l^RniT) = Ro{T) + -fRi{T) + -f^R2{T) + .... (42) 

n=0 



We note that the degree of the polynomials Qn{T) for n > 1 is at most iV — 1. For 7 = 0, we know that Eq — and 
that this maximal eigenvalue is obtained for yi — 0. Therefore, we have 

Qo{T) = T^ and Ro{T) = {1 ~ xT)^ + x'^ {I - T)^ . (43) 

By substituting the power series (|^T|) and (|42p in the QR-equation we obtain a hierarchical system of linear 
equations for the polynomials i?„(T) and Qn{T). This system can be solved order by order by using the known 'initial 
conditions' Qq{T) and RqIT). 

We now solve the QR-equation to the first and second orders. 



A. First order calculation 



At first order, the QR-equation becomes 

Qi(T) [(1 - xT)^ +x^{l~ T)^] + T^Ri{T) = (1 - T)^Qi{xT) + (1 - xT)^x^Qi (^^ ) + Lx^ {1 - T)^T^ , (44) 
and the auxiliary equation ([55)1 becomes 



It is simpler to define the polynomial 



l(1)-x^Qi(-) =-iV. (45) 



i?i(T) = gi(T)-a;^gi(-V (46) 



and to rewrite equations (|44|) and ((45|l as follows 

(1 - xT)^Bi{T) - (1 - TfBi{xT) = {Lx^ {I - T)^ - i?i(T)) . (47) 

Bi{l) = -N . (48) 

Because Bi(T) and Qi{T) are of degree < A'' — 1 and noting that the term on the r.h.s. of equation (|T7|) is divisible 
by T^, we can reduce this equation modulo and write 

(1 - xT)^Bi{T) - (1 - T)^Bi{xT) = [T^] . (49) 

This equation allows to determine the polynomial Bi{T) up to a multiplicative constant /3o 

N-l , . 

B,{T) EE /3o(l - T)^ [T^], z.e., i?i(T) = /Jq J] (-1)'^ U T ' ■ (^0) 
The constant /3o is fixed using equation (|48p . Using the binomial identity (jB3|) . we find 



fc=0 

From this relation it follows that 



L 

N-l 



Using this formula and equation (|47p the following exact expression for Ri{T) is obtained: 



Af-l L 



(53) 



All negative powers of T in the above expression cancel out for the following reason: the coefficient of a term of the 
type T"*^ with rf > is obtained by imposing the condition p + r — N — d to the double sum in equation ((55)) . Because 
of this condition, the indices p and r can vary only from to A^ — d and they both have the same effective range. The 
sum in equation ([53|l is antisymmetric with respect to p and r and therefore it vanishes. This proves that Ri{T) is 
indeed a polynomial. 



B. Second order calculation 



At second order, the polynomial B2{T) defined as 



B2{T)^Q2{T)-x''Q2{-] , (54) 



X 



satisfies the following equation 

(1 - xr)^B2(r) - (1 - T)'^B2{xT) = L{1 - T)'^Qi{xT) - R^{T)Q^{T) + ( ^{1 - T)^ - i?2(r) ) • (55) 



If we write this relation modulo we obtain the simpler equation 

(1 - xTfB^iT) - (1 - TfB^ixT) = L(l - TfQi{xT) - i?i(T)Qi(T) [T^] , (56) 

where the expressions for Qi{T) and Ri{T) are given in equations ((52|) and ([53|) respectively. At order 2, the auxiliary 
equation ((5^ becomes 

B2{1) = Q2(l) - x^Q2 = -A^Qi(l) - ^ . (57) 

The polynomial B2 {T) is the sum of a special solution B2 (T) of equation ([5^)1 and of a term that is proportional to 
Bi{T), the solution of the homogeneous equation ([H]), i.e., 

B2{T)^B2{T) + CBi{T). (58) 

The proportionality constant C is fixed by using the auxiliary equation ([57]) . which leads to 

C=^ + Qi(l) + f, (59) 



where we have used -Bi(l) = from equation 

A special solution to the polynomial equation ([56|) is given by 



N-l 



B2{T) = ^(l-a;^-^-)gf^T\ (60) 



k=0 



L \ L \ ^ ( L \ ( L 



with ^ = ^ ^ < y ^ + y ^ -WT—^ )■ ■ (61) 

N 

The main steps to derive equation ((^T|) are given in Appendix]^ Finally, the polynomial Q2iT) is given by the linear 
combination 

Af-l Af-1 

g2(T)=^gfT'^ + C^g[.^V\ (62) 

fc=0 k=0 

where the constant C is given in equation (|59p. 

VI. EXACT FORMULAE FOR THE MEAN CURRENT AND ITS FLUCTUATIONS 

Solving the Q-R equation allows us to calculate the expansion of the largest eigenvalue i?max(7), order by order, 
and to calculate the cumulants of the total current. The largest eigenvalue i?max(7), can be expanded with respect 
to the parameter 7 as follows 

00 

£^max(7) = P(l -X)Y, Y'Er, . (63) 
n=0 



Using equations (jiU)) and the expansion of i?max(7) is given by 

+ 7 7^ rT7--7^ IvL +7 -777^ ZXl + T, ZXl - T, ZXl + ■ ■ ■ (64) 



(65) 



p{l-x) 1-x (1-.t)^ '\{1-x)L {l-x)^J ' \ 2{l-x)^ {l-x)^ 
From the expression (|43p for Ro{T), we find that 

(l-x)^ 1-x 
and we verify that the zcroth-order term Eq in i?max(7) vanishes. 

A. Calculation of the Current 

The current J, defined in equation corresponds to the coefficient of 7 in the expansion of i?niax(7)- To determine 
R'i{l), we start with equation ()37|) and expand it to the first order in 7: 

^'^^^ =-N' + -^( Qi(l) - x^Qi (-]]- QUI) + ^^-^Q'l ( -) = -N' + -^i?i(l) - B[{1) , (66) 



{1 — x)^ l — x\ \^ J J \^ J 1 

where in the last equahty we have used the definition of Bi{T) as given by equation (jUJ. We know that i?i(l) = ~N 
from equation (|48|) : the value of B[{1) is readily obtained from the expression of -Bi(T) given in equations ((50|l 
and 

Bi(l) = -i^^- (67) 



Thus, we have 



1^ ^ ^ + • (68) 



(1-x)^ VI -a; i-1 

Substituting this expression in the coefficient of 7 in equation (j64p , we find that the total current is given by 

This value agrees, of course, with the known formula, which is obtained very simply by using the fact that all the 
stationary configurations of ASEP on a ring are equiprobable. We recall that J represents the total current in the 
system; the current through a bond is given by J / L. Using the Bethe Ansatz to find J is certainly a very complicated 
and distorted way to retrive a back-of-an-envelope calculation. However, J is one of the simplest quantity associated 
with ASEP and the fact that nobody could extract such an elementary formula from the Bethe equations has been a 
standing puzzle for a long time. 

B. Calculation of the Diffusion constant 

The second order term in the perturbative expansion ([M]) allows us to calculate the diffusion constant. Indeed, 
thanks to equation p^ . we find that A ~ 2p(l — x)E2. Therefore, we have 

A 9 n ^'-^o(l) ^ ^-^i(l) ,7nA 

Hence, in order to calculate A, we also need i?2(l), which is determined in Appendix IbI After gathering all relevant 
terms, we are finally lead to the exact formula for the diffusion constant of the total current for the partially asymmetric 
exclusion process on a ring: 



N 



A = 



2p{l-x)L ^ ^l + x"- f L \f L 



[L-l] 



N 



"^^^l-x^ [ N + r ) [ N-r ) ' 



This formula agrees, of course, with the one obtained using the Matrix Representation method 34] (in that work, 
the fluctuations of the current through a bond were calculated exactly i.e., A/L^). From this exact expression, it is 
possible to deduce by finite size scaling that a tagged particle in an infinite system exhibits an anomalous diffusive 
behaviour with exponent 1/3 (instead of one 1/2). By taking the continuous limit L — )■ cx3 of equation (j7ip in the 
weakly asymmetric regime x ^ 1, with scaling variable (f> = {1 — x)^/L, it is possible to derive a scaling function 
for the KPZ equation that describes the cross-over from the linear Edwards- Wilkinson regime to the non-linear KPZ 
regime. We refer for more details to [3^ . 

We ernphasize that the calculation of A with the Bethe Ansatz is of the same order of complexity as with the Matrix 
method [34| but it is much simpler mathematically. The Bethe Ansatz requires only elementary mathematical objects 
such as polynomials and involves systematic calculations, whereas for the Matrix Ansatz one has to find (guess) a 
suitable algebra, prove that this algebra solves the problem and then evaluate traces of various operators requiring 
the use of remarkable identities on q-binomials [s^l • 

Furthermore, to calculate the higher cumulants of the current, one has to solve the QR-equation (|33p to the suitable 
order in 7. By contrast, there is absolutely no clue on how to extend the Matrix method to calculate, for example, 
the third cumulant of the current: the form of the algebra involved (if such an algebra does exist) is totally unknown. 



VII. CONCLUSION 



Most of the analytical studies of the ASEP are based on two different techniques, the Matrix Product method and 
the Bethe Ansatz. The Matrix representation is suitable to calculate stationary state observables, such as correlations, 
phase diagrams etc... A major drawback of this method is that there is no constructive method to generate matrices 
that are suitable for a given stochastic model: one has to rely on educated guesses, after some trials and errors. 
Nevertheless, the Matrix method, when applicable, is efficient and allows to derive elegant combinatorial results for 
finite systems. On the contrary, the Bethe Ansatz is a systematic procedure with such a wide range of applicability 
that it has grown into a subfield of theoretical physics: the theory of integrable systems. There exists a priori 
conditions, such as the Yang-Baxter relation, that insure that a system is integrable {i.e., it can be analyzed by Bethe 
Ansatz). Many methods have been developed to cope with the Bethe equations [s^, [3§|. However, it is very difficult 
to extract information for finite systems from the Bethe equations and usually one has to analyze these equations in 
the thermodynamic limit. 

For the TASEP, the Bethe equations have a fundamental decoupling property that has lead to many exact results 
[m, [13, HE] and in particular to the calculation of an exact formula for the large deviation function [3l|, [12] ■ For the 
partially asymmetric case, the Bethe equations are strongly coupled and therefore they have been rarely used. The 
only exact results derived from them were obtained by Kim et al. in the limit of an infinite size [2^.[23|. In this paper, 
we have been able to overcome this technical difficulty thanks to a reformulation of the Bethe equations as a mere 
problem of polynomial divisibility that can be solved perturbatively in the fugacity parameter. We have calculated 
the mean value J of the current and its fluctuations A. Obviously, the calculation of J from Bethe Ansatz is much 
more difficult than the elementary derivation. However, the calculation of A with the Bethe Ansatz is less difficult 
than that with the Matrix method [13]. Furthermore, the perturbative analysis of the Bethe Ansatz can be extended 
a priori to any order to derive higher cumulants of the current. It is not known if the Matrix method can be applied 
to such calculations. 

The reformulation of the Bethe equations that we used here, is akin to the functional Bethe Ansatz [H,[H,[1§]. This 
method can be generalized to many other problems: higher moments of the ASEP current (S. Prolhac, in preparation), 
subleading correction to the large deviation function of the symmetric exclusion process, systems with different classes 
of particles. We also believe that the method followed here could be applied to the ASEP with open boundaries for 
which the Bethe equations have been derived recently (iol . [4ll ] . For the open TASEP with all rates equal to one, it 
is known from the Matrix method that the mean stationary current is given by the ratio of two consecutive Catalan 
numbers [13]: can this rather simple result be derived from Bethe Ansatz? 



Acknowledgments 



We thank Olivier Golinelli for many helpful discussions and S. Mallick for a careful reading of the manuscript. 



APPENDIX A: DERIVATION OF EQUATION ((6T|) 



We want to derive the formula ((M]) for B2{T) which is a particular solution of equation ((5^ . We substitute the 
formal expression (|60p of B2{T) in equation ([K^ and use the known explicit formulae for Qi{T) and Ri{T) (given 
in equations (|52p and (|53p . respectively). After identifying the terms of the same degree in T, the following linear 
system of equations is obtained: 

k+p=m ^ ^ k+p+r=m+N ^ / k+p=m \ / V / 

0<fe<Af-l 0<fc,p<Af-l 0<A:<Af-l 

0<p<L 0<r<L 0<p<L 



where we have introduced 



L 

c, = (-1)^+^+1 (A2) 

The system (|Aip is a triangular system of N equations, parametrized by the integer m, with 0<m<7V— l.In 
equation (|A1[) . we have written explicitely the ranges for all the dummy variables. However, some pieces of information 
are redundant: for example, we know that < m < A'^ — 1; therefore, if fc + p = m then both k and p must lie between 
and TV — 1 (recall that a binomial coefficient with a negative entry is equal to 0). In the following, we shall not write 
such superfluous information. 

We now start by transforming the r.h.s. of equation (jAip : we notice that the first sum on the r.h.s. is formally 
antisymmetric with respect to the indices r and p. However, this sum does not vanish identically because the range 
of these two variables is not the same. If the range of r were from to A^ — 1, the total sum would be equal to zero. 
In other words, the terms in the range 0<r<A — Ido not contribute to the sum, only the terms with N < r < L 
contribute. This sum is thus given by 

L \ f L \ f L\ xP -x' L \ V- f L \ f L 



r=N Q<k,p<N-l ^ / ^ / ^ / ,-=0 ^ ^ fe+p+r=m ^ ^ ^ 

k+p+r=m+N 0<k,p<N-l 



L-N 



^EIa + t) E ( p ) ( fc ^ r ) + 1 _ j;N+r^k ) (A3) 

where, we have first replaced the dummy variable r by r — A and then, to derive the last equality, we use the identity 
(xf — x^'^'^)/{l — X'^~^) = x^'^^ + [x^ — x^'^'^)/{l — x^"^) and replace A: by fc — r. Thus, we rewrite the r.h.s. of 



equation (|A1|) as follows: 



\ I \ \ I \ I \ I ^k , \ I \ \ I \ I \ ^k 



/^\N + r) E ( p ) ( fc-r ) 1^''+ l-x^+'-fe ) (a) E \ pl\k 
E(A + r) E ip)ik-r)^''^^[N + r) E ip)ik-r) I - a;W+'-fc ' (A'^) 



r— 1 ^ ^ k-\-p—m ^ ^ / r— ^ ^ k-\-p—m 

The first term in the last equality is now rewritten using the following identity 

k+p=m ^ / ^ / k+p=in ^ / ^ / 

[This identity is readily proved after multiplying both sides by (1 — a;''), cancelling the common x'^^^ term and noticing 
that the remaining terms are identical up to a notation change.] 
Finally, the initial system (jAip becomes: 



[x^ - x'')ck = (A6) 
EIa + t) E ( p ) ( fc- r ) + E ( A + r ) E [pjik- 



k-^-p—m ^ r— fc+p— m 



Clearly, the solution of this equation is given by 



L \ f L \ ( L \ I L 



= + 2^ — 1 _ .^N+r-k — • y^") 

r—l r— 

This ends the proof of the formula (|6ip . 

APPENDIX B: SOME USEFUL STEPS IN THE CALCULATION OF A 

1. Binomial formulae 

In the sequel, we shall use repeatedly the following elementary binomial formulae: 







L-l 
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^ L-l 
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iL-p)[ : ) = i( . ) • (B2) 
L \f L \ L^if L-l \ f L-l\ fL-l\f L-l 



N + r \ N -r 2'^\\N + r-l \N-r \N + r \N-r-l 



L f L-l \ f L-l \ _ N{L -N) ( L 



2 \ N -I \ N 2L \ N 
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(B4) 



2. An expression for the diffusion constant 



We start with equation ([57)1 and expand the polynomials Q and R to the second order in 7. This allows us to derive 
the following expression for i?2(l)- We obtain 

go(l)i?i(i) + gi(i)i?'i(l) + Q2{l)R'oil) + Qo(l)^2(l) + Q'i(l)i?i(l) + Q2(l)i?o(l) 

= -Lx^'+'il ~ x)^-ig2 (^i) + - xfQ', (^i^ . (B5) 

We know from equation (|43|) that Qo(l) = 1, Qoi^) = ^o(l) = (1 ~ ^)^- From equation ([40|). we deduce 
= N{1 - x)^ and i?2(l) = iV^(l - a;)^/2. Finally, equations ^ and dM]) give the values of and R'^{1). 

We also use equation ((57)) to express Q2(l/a;) in terms of (5i(l) and Q2(l)- Substituting this information into the 
previous expression leads to (remark that terms proportional to Q2(l) cancel out): 



iV^^Qi(l) - NQ'Al) - g^(l) + x^'-'Q'^ . (B6) 



{l~x)^ 2 2(l-x) ' " L 

Inserting this expression into the formula (jTU]) for A gives 

. iX! _ _ iv|^0,(i) + iVOKD + Bid) , (B7) 

where we have used the definition ((54|) of the polynomial B-ziT). With the help of equations (f58|) ([59]) . and ([67]) we 
get 

5^(1) = i?2'(l) + B[{1) + Qi(l) + I j = B2'(l) - i^B2(l) - LiV^Qi(l) - i^'^II^T) ■ 



Substituting this expression in equation (|B9p . we conclude that 

The values of all the terms that appear in this equation are known. We now evaluate each of these terms separately. 

3. Calculation of some exact expressions 

The value of Qi{l) is easily obtained from the expression ([5^ of Qi{T): 

Similarly, we have 

Q'AD-rh^ ^ (BID 

To calculate i?2(l), we start from the formula for B2{T) given in equations ([50]) and (pT|) 




(B12) 



Exchanging the double sum and using equations (jBip and (|B3p . we rewrite the first term on the r.h.s of this expression 
as follows: 



k=0 



L 

N -r 



In the second term on the r.h.s. of (|B12p . we remark that the effective range of the variable r is from to fc and 
we replace r by r' = fc — r. We then transform this term in a manner similar to that described in equation (|B13p . 
Finally, the expression (|B12p simplifies to: 

L ^x'-{N-r)+N + r f L \ f L \ LN ^""^^'^ ( TV - r 



Using similar steps, we find that B2 (1) is given by 

L ^x'-{N-r){LN^r-L) + {N + r){LN + r-L) f L \ f L 

^2(1) = ryr^Z^ 



l-x'' \ N + r J \ N - r 



N 



LN Ni-mLN-Lr + r-L) 



(L-l) 



L 



N 



To conclude our calculation, we must substitute equations (jBlOl) (IB11[) . (|B14p . and (jBlSp into the formula (|B9p for 
the diffusion constant. We find that all the terms that contain only one binomial factor i.e., terms proportional to 

(—1)'' ^ ^ ^ /(f — x'^) cancel out amongst themselves. After some elementary simplications, we are left with 
A 2L nl + x"^ / L \ / L \ 2LN f L \f L \ 



(B16) 

The last two terms cancel with each other according to the identity (|B4p . This ends the proof of equation (|7T|) . 



APPENDIX C: FUNCTIONAL BETHE ANSATZ FOR TASEP 

We consider here the special case of the TASEP (which corresponds to p = 1 and q — x ~ {)). We exp lain how 
to retrieve from the Q-R equation p3p the parametric representation of £'max(7) that was obtained in |3l| by using 
contour integrals. For a; = 0, the functional equation reduces to 



Q{T)R{T)^T^ + [-l)^-^B{l-T)^ with B = [-l)^-^e^''Q[Q) . (CI) 
From equation (j43p . we find that the zeroth order polynomials for the TASEP are simply given by 

Qo(r) = r^ and i?o(T) = l. (C2) 

The perturbative expansions ((1T|) and (jH]) can be rewritten as 

Q(r) = r^ + 7Q(r) and i?(r) = 1 + 77^(^) , (ca) 

where Q{T) is a polynomial of degree N — 1 and TZ{T) is of degree L ~ N (the coefficients of these two polynomials 
are functions of 7). We note, in particular, that Q(0) is of order 7 and thus B defined in equation (|C1|) is also of 
order 7 and is a small parameter. Dividing both sides of equation (jCip by and taking the logarithm, we obtain 



^ ^ ^ k=l 

where we have developped the logarithm in powers of B. We remark that the r.h.s. of this equation is a series that 
contains both positive and negative powers of T. But, equation (|C3p implies that Q{T)/T^ = 1 + ^Q{T) /T^ , i.e., 
Q{T)/T^ is a polynomial in the variable 1/T of degree N. Therefore, the expansion of log {Q{T) /T^^ w.r.t. 7 (or 
B) can only generate negative powers of T. Similarly, from equation (|C3p we have log_R(r) = log(l + ^Ti{T)) and 
the expansion of this term can generate only positive powers of T. Therefore, the identification between the l.h.s. and 
the r.h.s. of equation (|C4p is unique and we have 

.o.(^)^|^^i^|Vir(-)T-~ ,C5) 

logi?(T) = k ^ (-1)' ( 7 ) 

fe=l j=kN ^ ^ 

For the TASEP, equations jMl) and ^ reduce to 

i?max(7) = -f|^--^logi?(r)|^^^ and 7=llogi?(l). (C7) 
From equation (IC6p . we obtain (with the help of equations (IBip and (jB3p to calculate the sums over j): 



{kNy.ikL - kNy. ' ^ ' 



These two equations are precisely those derived in [3l|. They provide a parametric formula for -Ernax(7) that allows 
to calculate the large deviation function of the current and its cumulants to any required order. 
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